﻿using System;

namespace Atomic.Libraries.Mathematics
{
	public static partial class SpecialFunctions
	{
		/// <summary>
		/// The polylogarithm $Li_n(x)$. Real part of the polylogarithm for x > 1.
		/// </summary>
		public static double PolyLog(int n, double x)
		{
			// Implementation using Cephes public domain library:
			// http://www.netlib.org/cephes/
			//
			// Translated from old C-code:
			// http://www.netlib.org/cephes/misc.tgz
			// http://www.netlib.org/cephes/cmath.tgz

			if (x > 1.0)
			{
				// Use real part of the polylogarithm here.
				return Complex.Re(PolyLog(n, (Complex)x));
			}

			return PolyLogImplementation.polylog(n, x);
		}

		/// <summary>
		/// The general complex-valued polylogarithm $Li_n(z)$.
		/// </summary>
		public static Complex PolyLog(int n, Complex z)
		{
			return PolyLogImplementation.cpolylog(n, z);
		}

		private static class PolyLogImplementation
		{
			private static double machep, pi;
			private static double[] polylog_a4, polylog_b4, spence_a, spence_b, zetac_azetac, zetac_p, zetac_q, zetac_r, zetac_s, zetac_a, zetac_b;

			static PolyLogImplementation()
			{
				machep = 1.11022302462516e-16;
				pi = Math.PI;
				polylog_a4 = new double[] { 0.030561449220894907, 0.32430864841625817, 0.28778472814618761, 0.070912677858861808, 0.0064664600724566208, 0.00024502330192965427, 4.0316553646277049e-06, 2.8841691639094679e-08, 8.6800670024665943e-11, 1.025983405866371e-13, 4.2334683135382725e-17, 4.9594220350662069e-21, 1.0593658675852757e-25 };
				polylog_b4 = new double[] { 2.8212624036003109, 1.7802211248813271, 0.37788882118678757, 0.031938870400743377, 0.0011612524184980964, 1.8673623748298707e-05, 1.3190227797152944e-07, 3.9427552565556031e-10, 4.6443269689863969e-13, 1.913336021014307e-16, 2.2400418146260698e-20, 4.784036597230791e-25, 0.0 };
				spence_a = new double[] { 4.6512858607399003e-05, 0.0073158904523809472, 0.13384763957830903, 0.8796913117545303, 2.7114985119655346, 4.2569715600812179, 3.2977134098522511, 1.0 };
				spence_b = new double[] { 0.00069099048891255329, 0.025404376393254439, 0.2829748606025681, 1.4117259775183106, 3.6380053334513707, 5.0327888014331696, 3.5477134098522511, 1.0 };
				zetac_azetac = new double[] { -1.5, 1.7976931348623157e+308, 0.64493406684822641, 0.20205690315959429, 0.082323233711138186, 0.036927755143369927, 0.01734306198444914, 0.0083492773819228271, 0.0040773561979443396, 0.0020083928260822143, 0.00099457512781808526, 0.00049418860411946453, 0.00024608655330804832, 0.00012271334757848915, 6.1248135058704828e-05, 3.0588236307020493e-05, 1.5282259408651871e-05, 7.6371976378997626e-06, 3.8172932649998402e-06, 1.908212716553939e-06, 9.5396203387279621e-07, 4.7693298678780645e-07, 2.38450502727733e-07, 1.1921992596531106e-07, 5.960818905125948e-08, 2.9803503514652279e-08, 1.4901554828365043e-08, 7.4507117898354301e-09, 3.7253340247884573e-09, 1.8626597235130491e-09, 9.3132743241966817e-10 };
				zetac_p = new double[] { 5.85746514569725319540e11, 2.57534127756102572888e11, 4.87781159567948256438e10, 5.15399538023885770696e9, 3.41646073514754094281e8, 1.60837006880656492731e7, 5.92785467342109522998e5, 1.51129169964938823117e4, 2.01822444485997955865e2 };
				zetac_q = new double[] { 3.90497676373371157516e11, 5.22858235368272161797e10, 5.64451517271280543351e9, 3.39006746015350418834e8, 1.79410371500126453702e7, 5.66666825131384797029e5, 1.60382976810944131506e4, 1.96436237223387314144e2 };
				zetac_r = new double[] { -3.28717474506562731748e-1, 1.55162528742623950834e1, -2.48762831680821954401e2, 1.01050368053237678329e3, 1.26726061410235149405e4, -1.11578094770515181334e5 };
				zetac_s = new double[] { 1.95107674914060531512e1, 3.17710311750646984099e2, 3.03835500874445748734e3, 2.03665876435770579345e4, 7.43853965136767874343e4 };
				zetac_a = new double[] { 8.70728567484590192539e6, 1.76506865670346462757e8, 2.60889506707483264896e10, 5.29806374009894791647e11, 2.26888156119238241487e13, 3.31884402932705083599e14, 5.13778997975868230192e15, -1.98123688133907171455e15, -9.92763810039983572356e16, 7.82905376180870586444e16, 9.26786275768927717187e16 };
				zetac_b = new double[] { -7.92625410563741062861e6, -1.60529969932920229676e8, -2.37669260975543221788e10, -4.80319584350455169857e11, -2.07820961754173320170e13, -2.96075404507272223680e14, -4.86299103694609136686e15, 5.34589509675789930199e15, 5.71464111092297631292e16, -1.79915597658676556828e16 };
			}

			public static double polylog(int n, double x)
			{
				double h, k, p, s, t, u, xc, z;
				int i, j;

				// This recurrence provides formulas for n < 2:
				//
				//  d                 1
				//  --   Li (x)  =   ---  Li   (x)  .
				//  dx     n          x     n-1

				if (n < -1)
				{
					// Trivial expressions for n < -1 isn't implemented.
					throw new NotImplementedException();
				}

				if (n == -1)
				{
					p = 1.0 - x;
					u = x / p;
					s = u * u + u;
					return s;
				}

				if (n == 0)
				{
					s = x / (1.0 - x);
					return s;
				}

				// Not defined for x > 1. Use cpolylog if you need that.
				if (x > 1.0)
				{
					return double.NaN;
				}

				if (n == 1)
				{
					s = -log(1.0 - x);
					return s;
				}

				// Argument +1
				if (x == 1.0 && n > 1)
				{
					s = zetac(n) + 1.0;
					return s;
				}

				// Argument -1:
				//                    1-n
				// Li (-z)  = - (1 - 2   ) Li (z)
				//   n                       n
				if (x == -1.0 && n > 1)
				{
					// Li_n(1) = zeta(n)
					s = zetac(n) + 1.0;
					s = s * (powi(2.0, 1 - n) - 1.0);
					return s;
				}

				// Inversion formula:
				//                                                  [n/2]   n-2r
				//               n                  1     n           -  log    (z)
				// Li (-z) + (-1)  Li (-1/z)  =  - --- log (z)  +  2  >  ----------- Li  (-1)
				//   n               n              n!                -   (n - 2r)!    2r
				//                                                   r=1
				if (x < -1.0 && n > 1)
				{
					double q, w;
					int r;

					w = log(-x);
					s = 0.0;
					for (r = 1; r <= n / 2; r++)
					{
						j = 2 * r;
						p = polylog(j, -1.0);
						j = n - j;
						if (j == 0)
						{
							s = s + p;
							break;
						}
						q = (double)j;
						q = pow(w, q) * p / fac(j);
						s = s + q;
					}
					s = 2.0 * s;
					q = polylog(n, 1.0 / x);
					if ((n & 1) == 1)
					{
						q = -q;
					}
					s = s - q;
					s = s - pow(w, (double)n) / fac(n);
					return s;
				}

				if (n == 2)
				{
					if (x < 0.0 || x > 1.0)
					{
						return spence(1.0 - x);
					}
				}

				// The power series converges slowly when x is near 1.  For n = 3, this
				// identity helps:
				//
				// Li (-x/(1-x)) + Li (1-x) + Li (x)
				//   3               3          3
				//                2                               2                 3
				//  = Li (1) + (pi /6) log(1-x) - (1/2) log(x) log (1-x) + (1/6) log (1-x)
				//      3
				if (n == 3)
				{
					p = x * x * x;
					if (x > 0.8)
					{
						// Thanks to Oscar van Vlijmen for detecting an error here.
						u = log(x);
						s = u * u * u / 6.0;
						xc = 1.0 - x;
						s = s - 0.5 * u * u * log(xc);
						s = s + pi * pi * u / 6.0;
						s = s - polylog(3, -xc / x);
						s = s - polylog(3, xc);
						s = s + zetac(3.0);
						s = s + 1.0;
						return s;
					}
					// Power series
					t = p / 27.0;
					t = t + 0.125 * x * x;
					t = t + x;

					s = 0.0;
					k = 4.0;
					do
					{
						p = p * x;
						h = p / (k * k * k);
						s = s + h;
						k += 1.0;
					} while (fabs(h / s) > 1.1e-16);
					return s + t;
				}

				if (n == 4)
				{
					if (x >= 0.875)
					{
						u = 1.0 - x;
						s = polevl(u, polylog_a4, 12) / p1evl(u, polylog_b4, 12);
						s = s * u * u - 1.202056903159594285400 * u;
						s += 1.0823232337111381915160;
						return s;
					}
					goto pseries;
				}

				if (x < 0.75)
				{
					goto pseries;
				}

				// This expansion in powers of log(x) is especially useful when
				// x is near 1.
				//
				// See also the pari gp calculator.
				//
				//                   inf                  j
				//                    -    z(n-j) (log(x))
				// polylog(n,x)  =    >   -----------------
				//                    -           j!
				//                   j=0
				//
				//   where
				//
				//   z(j) = Riemann zeta function (j), j != 1
				//
				//                           n-1
				//                            -
				//   z(1) =  -log(-log(x)) +  >  1/k
				//                            -
				//                           k=1

				z = log(x);
				h = -log(-z);
				for (i = 1; i < n; i++)
				{
					h = h + 1.0 / i;
				}
				p = 1.0;
				s = zetac(n) + 1.0;
				for (j = 1; j <= n + 1; j++)
				{
					p = p * z / j;
					if (j == n - 1)
					{
						s = s + h * p;
					}
					else
					{
						s = s + (zetac(n - j) + 1.0) * p;
					}
				}
				j = n + 3;
				z = z * z;
				while (true)
				{
					p = p * z / ((j - 1) * j);
					h = (zetac(n - j) + 1.0);
					h = h * p;
					s = s + h;
					if (fabs(h / s) < machep)
					{
						break;
					}
					j += 2;
				}
				return s;

			pseries:
				p = x * x * x;
				k = 3.0;
				s = 0.0;
				do
				{
					p = p * x;
					k += 1.0;
					h = p / powi(k, n);
					s = s + h;
				} while (fabs(h / s) > machep);
				s += x * x * x / powi(3.0, n);
				s += x * x / powi(2.0, n);
				s += x;
				return s;
			}

			private static double spence(double x)
			{
				double w, y, z;
				int flag;

				if (x < 0.0)
				{
					return double.NaN;
				}

				if (x == 1.0)
				{
					return 0.0;
				}

				if (x == 0.0)
				{
					return pi * pi / 6.0;
				}

				flag = 0;

				if (x > 2.0)
				{
					x = 1.0 / x;
					flag |= 2;
				}

				if (x > 1.5)
				{
					w = (1.0 / x) - 1.0;
					flag |= 2;
				}
				else if (x < 0.5)
				{
					w = -x;
					flag |= 1;
				}
				else
				{
					w = x - 1.0;
				}

				y = -w * polevl(w, spence_a, 7) / polevl(w, spence_b, 7);

				if ((flag & 1) == 1)
				{
					y = (pi * pi) / 6.0 - log(x) * log(1.0 - x) - y;
				}

				if ((flag & 2) == 2)
				{
					z = log(x);
					y = -0.5 * z * z - y;
				}

				return y;
			}

			private static double polevl(double x, double[] coef, int n)
			{
				double ans = coef[0];
				for (int i = 1; i <= n; i++)
				{
					ans = ans * x + coef[i];
				}
				return ans;
			}

			private static double p1evl(double x, double[] coef, int n)
			{
				double ans = x + coef[0];
				for (int i = 1; i < n; i++)
				{
					ans = ans * x + coef[i];
				}
				return ans;
			}

			private static double pow(double x, double y)
			{
				return Math.Pow(x, y);
			}

			private static double powi(double x, int n)
			{
				return Math.Pow(x, n);
			}

			private static double log(double x)
			{
				return Math.Log(x);
			}

			private static double fac(int i)
			{
				return SpecialFunctions.Gamma(i + 1.0); // FIXME
			}

			private static double fabs(double x)
			{
				return Math.Abs(x);
			}

			private static double zetac(double x)
			{
				int i;
				double a, b, s, w;

				if (x < 0.0)
				{
					if (x < -170.6243)
					{
						return double.NaN;
					}
					s = 1.0 - x;
					w = zetac(s);
					b = sin(0.5 * pi * x) * pow(2.0 * pi, x) * gamma(s) * (1.0 + w) / pi;
					return b - 1.0;
				}

				if (x >= 127)
				{
					return 0.0; // because first term is 2^-x
				}

				// Tabulated values for integer argument
				w = floor(x);
				if (w == x)
				{
					i = (int)x;
					if (i < 31)
					{
						return zetac_azetac[i];
					}
				}

				if (x < 1.0)
				{
					w = 1.0 - x;
					a = polevl(x, zetac_r, 5) / (w * p1evl(x, zetac_s, 5));
					return a;
				}

				if (x == 1.0)
				{
					return double.PositiveInfinity;
				}

				if (x <= 10.0)
				{
					b = pow(2.0, x) * (x - 1.0);
					w = 1.0 / x;
					s = (x * polevl(w, zetac_p, 8)) / (b * p1evl(w, zetac_q, 8));
					return s;
				}

				if (x <= 50.0)
				{
					b = pow(2.0, -x);
					w = polevl(x, zetac_a, 10) / p1evl(x, zetac_b, 10);
					w = exp(w) + b;
					return w;
				}

				// Basic sum of inverse powers

				s = 0.0;
				a = 1.0;
				do
				{
					a += 2.0;
					b = pow(a, -x);
					s += b;
				} while (b / s > machep);

				b = pow(2.0, -x);
				s = (s + b) / (1.0 - b);
				return s;
			}

			private static double gamma(double x)
			{
				return SpecialFunctions.Gamma(x);
			}

			private static double exp(double x)
			{
				return Math.Exp(x);
			}

			private static double floor(double x)
			{
				return Math.Floor(x);
			}

			private static double sin(double x)
			{
				return Math.Sin(x);
			}

			public static Complex cpolylog(int n, Complex x)
			{
				Complex h, p, s, u, z;
				double k;
				double ah, as_;
				int i, j;

				// This recurrence provides formulas for n < 2.
				//
				// d                 1
				// --   Li (x)  =   ---  Li   (x)  .
				// dx     n          x     n-1

				if (n < -1)
				{
					// Trivial expressions for n < -1 isn't implemented.
					throw new NotImplementedException();
				}

				if (n == -1)
				{
					p = 1.0 - x;
					u = x / p;
					s = u * u + u;
					return s;
				}

				if (n == 0)
				{
					s = x / (1.0 - x);
					return s;
				}

				if (n == 1)
				{
					s = -clog(1.0 - x);
					return s;
				}

				// Argument +1
				if (x == (1.0 + 0.0 * Complex.I) && n > 1)
				{
					s = zetac(n) + 1.0 + 0.0 * Complex.I;
					return s;
				}

				// Argument -1.
				//                    1-n
				// Li (-z)  = - (1 - 2   ) Li (z)
				//   n                       n
				if (x == (-1.0 + 0.0 * Complex.I) && n > 1)
				{
					// Li_n(1) = zeta(n)
					s = zetac(n) + 1.0 + 0.0 * Complex.I;
					s = s * (powi(2.0, 1 - n) - 1.0);
					return s;
				}

				// Use inversion formula for large x.
				ah = fabs(creal(x)) + fabs(cimag(x));
				if (ah > 3.0)
				{
					return (invformula(n, x));
				}

				if ((fabs(creal(x)) + fabs(cimag(x))) < 0.75)
				{
					goto pseries;
				}

				// The expansion in powers of log(x) is especially useful when
				// x is near 1.
				//
				// See also the pari gp calculator.
				//
				//                   inf                  j
				//                    -    z(n-j) (log(x))
				// polylog(n,x)  =    >   -----------------
				//                    -           j!
				//                   j=0
				//
				// where
				//
				// z(j) = Riemann zeta function (j), j != 1
				//
				//                         n-1
				//                          -
				// z(1) =  -log(-log(x)) +  >  1/k
				//                          -
				//                         k=1

				z = clog(x);
				h = -clog(-z);
				for (i = 1; i < n; i++)
				{
					h = h + 1.0 / i;
				}
				p = 1.0;
				s = zetac(n) + 1.0;
				for (j = 1; j <= n + 1; j++)
				{
					p = p * z / j;
					if (j == n - 1)
					{
						s = s + h * p;
					}
					else
					{
						s = s + (zetac(n - j) + 1.0) * p;
					}
				}
				j = n + 3;
				z = z * z;
				while (true)
				{
					p = p * z / ((j - 1) * j);
					h = (zetac(n - j) + 1.0);
					h = h * p;
					s = s + h;
					ah = fabs(creal(h)) + fabs(cimag(h));
					as_ = fabs(creal(s)) + fabs(cimag(s));
					if (ah / as_ < machep)
					{
						break;
					}
					j += 2;
				}
				return s;

			pseries:
				s = 0.0;
				ah = fabs(creal(x)) + fabs(cimag(x));
				if (ah < 1.0e-6)
				{
					goto smseries;
				}

				p = x * x * x;
				k = 3.0;
				do
				{
					p = p * x;
					k += 1.0;
					h = p / powi(k, n);
					s = s + h;
					ah = fabs(creal(h)) + fabs(cimag(h));
					as_ = fabs(creal(s)) + fabs(cimag(s));
				} while (ah / as_ > 1.1e-16);

			smseries:
				s += x * x * x / powi(3.0, n);
				s += x * x / powi(2.0, n);
				s += x;
				return s;
			}

			private static Complex invformula(int n, Complex x)
			{
				// Inversion formula:
				//                                                  [n/2]   n-2r
				//               n                  1     n           -  log    (z)
				// Li (-z) + (-1)  Li (-1/z)  =  - --- log (z)  +  2  >  ----------- Li  (-1)
				//   n               n              n!                -   (n - 2r)!    2r
				//                                                   r=1
				Complex p, q, s, w, m1;
				int r, j;

				m1 = -1.0 + 0.0 * Complex.I;
				w = clog(-x);
				s = 0.0 + 0.0 * Complex.I;
				for (r = 1; r <= n / 2; r++)
				{
					j = 2 * r;
					p = cpolylog(j, m1);
					j = n - j;
					if (j == 0)
					{
						s = s + p;
						break;
					}
					q = (double)j + 0.0 * Complex.I;
					q = cpow(w, q) * (creal(p) / fac(j));
					s = s + q;
				}
				s = 2.0 * s;

				q = cpolylog(n, 1.0 / x);
				if ((n & 1) == 1)
				{
					q = -q;
				}
				s = s - q;

				s = s - cpow(w, n) / fac(n);
				return s;
			}

			private static double creal(Complex x)
			{
				return Complex.Re(x);
			}

			private static double cimag(Complex x)
			{
				return Complex.Im(x);
			}

			private static Complex clog(Complex x)
			{
				return Complex.Log(x);
			}

			private static Complex cpow(Complex x, Complex y)
			{
				return Complex.Exp(Complex.Log(x) * y); // FIXME
			}
		}
	}
}
